day 7: multilevel models

multilevel tadpoles

review (or introduction?)

(Ignore probability and estimation for now)

\[ Y_i = b_0 + b_1X_{1i} + ... + \epsilon_i \] This assumes each observation \((i)\) is independent of all other observations. But what if our data were clustered in some way? We might be able to express each observation \(i\) as belonging to a specific cluster, \(j\).

\[ Y_{ij} = b_0 + b_1X_{1ij} + ... + \epsilon_{ij} \]

But what is the problem with this?

multilevel models

We’re starting our unit on multilevel models, which can be thought of as models that “remember” features of clusters of data as they learn about all the clusters. The model will pool information across clusters (e.g., our estimates about cluster A will be informed in part by clusters B, C, and D). This tends to improve estimates about each cluster. Here are some other benefits of multilevel modeling:

  1. improved estimates for repeated sampling. If you try to fit a single-level model to these data, you’ll over- or under-fit the data.
  2. improved estimates for imbalance in sampling. prevent over-sampled clusters from dominating inference, while also balancing the fact that larger clusters have more information.
  3. estimates of variation. model variation explicitly!
  4. avoid averaging, retain variation. averaging manufactures false confidence (artificially inflates precision) and introduces arbitrary data transformations.

Multilevel modeling should be your default approach.

(Still non-Bayesian)

Intercept-only model

\[\begin{align*} \text{Level 1}&\\ Y_{ij} &= \beta_{0j} + \epsilon_{ij}\\ \text{Level 2}&\\ \beta_{0j} &= \gamma_{00} + U_{0j} \\ U_{0j} &\sim \text{Normal}(0, \tau^2_{00}) \\ \epsilon_{ij} &\sim \text{Normal}(0, \sigma^2) \\ \end{align*}\]

We can rewrite this as the “combined” model:

\[ Y_{ij} = \gamma_{00} + U_{0j} + \epsilon_{ij}\\ \] Level 1 is where you have data that repeats within your grouping or clustering data. Is your cluster classrooms? Then students are Level 1. Is your cluster people? Then observations are Level 1.

Level 2 takes the parameters at Level 1 and decomposes them into a fixed component \((\gamma)\) that reflects the average and, if desired, the individual deviations around that fixed effect \((U)\).

\[\begin{align*} \text{Level 1}&\\ Y_{ij} &= \beta_{0j} + \epsilon_{ij}\\ \text{Level 2}&\\ \beta_{0j} &= \gamma_{00} + U_{0j} \\ U_{0j} &\sim \text{Normal}(0, \tau^2_{00}) \\ \epsilon_{ij} &\sim \text{Normal}(0, \sigma^2) \\ \end{align*}\]

Each person’s score is therefore divided into the “fixed” and “random” components:

\[\begin{align*} Y_{i,16} &= \beta_{0, 16} + \epsilon_{ij} \\ ... \\ B_{0, 16} &= \gamma_{00} + U_{0,16} \end{align*}\]

(Random noise is assumed to vary constantly across observations across participants (in this model).)

example: multilevel people

iFit study

  • 193 participants recruited from a commercial weight loss study as well as the general population in two California counties
  • Wore a Fitbit Charge for 24 hours/day for at least 100 days
  • Daily diary survey about several variables, including affect (PANAS) resulting in a positive and negatve affect score for each day.
data_path = "https://raw.githubusercontent.com/sjweston/uobayes/refs/heads/main/files/data/external_data/williams.csv"
d <- read.csv(data_path)
rethinking::precis(d)
                 mean         sd       5.5%      94.5%      histogram
id        98.61029694 63.7493291 10.0000000 207.000000   ▇▇▇▃▅▅▅▃▃▃▂▂
female     0.57016803  0.4950710  0.0000000   1.000000     ▅▁▁▁▁▁▁▁▁▇
PA.std     0.01438236  1.0241384 -1.6812971   1.751466       ▁▁▃▇▇▃▁▁
day       44.17962096 27.6985612  6.0000000  92.000000     ▇▇▇▇▅▅▅▃▃▃
PA_lag     0.01992587  1.0183833 -1.7204351   1.717036     ▁▂▃▅▇▇▅▃▂▁
NA_lag    -0.04575229  0.9833161 -0.8750718   1.990468 ▇▃▂▁▁▁▁▁▁▁▁▁▁▁
steps.pm   0.05424387  0.6298941 -1.0258068   1.011356       ▁▂▇▇▃▁▁▁
steps.pmd  0.02839160  0.7575395 -1.1235951   1.229974 ▁▁▇▇▁▁▁▁▁▁▁▁▁▁
NA.std    -0.07545093  0.9495660 -1.0484293   1.811061     ▁▁▁▇▂▁▁▁▁▁
d %>% count(id) %>% rethinking::precis()
       mean       sd  5.5%  94.5%    histogram
id 120.6425 71.29547 11.56 227.44 ▇▇▇▅▇▇▇▇▇▇▇▇
n   67.5285 28.01953 17.00  99.00   ▁▁▂▂▁▂▂▂▂▇
Code 
set.seed(114)
sample_id = sample(unique(d$id), replace=F, size=6)
d %>%
  filter(id %in% sample_id) %>% 
  mutate(id = factor(id)) %>% 
  ggplot(aes(x = day, y = PA.std, group = id, fill = id)) + 
  geom_point(aes(color = id)) + 
  geom_line( aes(color = id), alpha=.5) +
  labs(y="Positive Affect") + 
  facet_wrap(~id) +
  theme(legend.position = "none")

What if we wanted to estimate each person’s conscientiousness score? One method would be to simply average scores for each person, but we lose a lot of information that way. Another option would be to treat each person as a group and model scores as a function of group. We can do this using an unpooled model.

\[\begin{align*} \text{PA.std}_i &\sim \text{Normal}(\mu_i,\sigma) \\ \mu_i &= \alpha_{\text{id}[i]} \\ \alpha_j &\sim \text{Normal}(0, 1.5) \text{ for }j=1,...,91 \\ \sigma &\sim \text{Exponential}(1) \end{align*}\]

d$id.f = as.factor(d$id)

m1 <- brm(
  data=d,
  family=gaussian,
  bf(PA.std ~ 0 + a,
     a ~ 0 + id.f,
     nl = TRUE),
  prior = c( prior(normal(0, 1.5), class=b, nlpar=a),
             prior(exponential(1), class=sigma)),
  iter=2000, warmup=1000, chains=4, cores=4, seed=9,
  file=here("files/models/m71.1")
)
m1
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: PA.std ~ 0 + a 
         a ~ 0 + id.f
   Data: d (Number of observations: 13033) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
a_id.f1       0.08      0.06    -0.04     0.21 1.00    11270     2736
a_id.f2       1.94      0.06     1.82     2.06 1.00    10974     2899
a_id.f3      -0.03      0.06    -0.15     0.09 1.00     8695     3022
a_id.f4      -0.31      0.06    -0.43    -0.20 1.00     8990     2238
a_id.f5       1.12      0.14     0.84     1.40 1.00    12665     2989
a_id.f6      -0.53      0.07    -0.66    -0.39 1.00     7837     2500
a_id.f7       0.19      0.07     0.05     0.32 1.00    10511     2437
a_id.f8       2.13      0.06     2.00     2.25 1.00    10646     2581
a_id.f9       1.04      0.08     0.89     1.20 1.00     9015     2893
a_id.f10      0.51      0.06     0.39     0.62 1.00     9612     2709
a_id.f11      0.09      0.06    -0.03     0.21 1.00     9234     2953
a_id.f12      0.97      0.06     0.85     1.09 1.00     8591     2734
a_id.f13      0.35      0.06     0.23     0.47 1.00    10446     3255
a_id.f14     -0.12      0.06    -0.24     0.01 1.00    10432     2906
a_id.f15     -0.48      0.06    -0.60    -0.36 1.00    10630     2845
a_id.f16      0.79      0.07     0.67     0.92 1.00     9938     2910
a_id.f17     -0.07      0.07    -0.19     0.06 1.00     9661     2624
a_id.f20     -0.60      0.06    -0.73    -0.48 1.00     9270     2729
a_id.f23     -0.58      0.06    -0.70    -0.46 1.00    10089     3045
a_id.f24      1.97      0.06     1.86     2.09 1.00     7588     2825
a_id.f25      0.29      0.06     0.17     0.41 1.00     8913     2177
a_id.f26     -0.43      0.06    -0.55    -0.31 1.00     9631     2987
a_id.f27      2.02      0.06     1.89     2.14 1.00     8945     2670
a_id.f28     -1.83      0.07    -1.96    -1.71 1.00     9091     2541
a_id.f29     -0.09      0.06    -0.22     0.03 1.00    12304     2658
a_id.f30      0.92      0.06     0.80     1.04 1.00    11099     2580
a_id.f31     -0.48      0.06    -0.60    -0.36 1.00    10211     2559
a_id.f32      1.34      0.06     1.21     1.46 1.00     9726     2749
a_id.f33     -0.78      0.06    -0.90    -0.66 1.00    10981     2697
a_id.f34      0.30      0.06     0.19     0.42 1.00     8492     2520
a_id.f35     -0.15      0.06    -0.27    -0.02 1.00    10766     2704
a_id.f36      0.09      0.06    -0.03     0.21 1.00    13422     2339
a_id.f37      0.73      0.06     0.61     0.85 1.00    11831     2532
a_id.f38     -1.52      0.06    -1.64    -1.41 1.00     8925     2423
a_id.f39     -0.15      0.06    -0.27    -0.04 1.00     8167     2733
a_id.f40      0.25      0.06     0.12     0.37 1.00    10679     2932
a_id.f41      0.78      0.06     0.65     0.90 1.00     8825     2950
a_id.f42     -0.45      0.06    -0.58    -0.33 1.00     9401     2473
a_id.f43      0.43      0.08     0.28     0.58 1.00     9432     2655
a_id.f44     -0.52      0.06    -0.64    -0.39 1.00    10960     2713
a_id.f46      0.18      0.07     0.05     0.31 1.00    10528     2871
a_id.f47     -0.52      0.06    -0.64    -0.41 1.00    10959     2772
a_id.f48     -0.22      0.07    -0.35    -0.10 1.00     7910     2716
a_id.f49     -1.36      0.06    -1.48    -1.24 1.00     8662     2579
a_id.f50     -0.79      0.06    -0.91    -0.68 1.00     9687     2657
a_id.f51      1.77      0.06     1.65     1.90 1.00     8890     2653
a_id.f52     -0.11      0.06    -0.23     0.01 1.00    10873     2560
a_id.f53      0.31      0.06     0.19     0.42 1.00     9033     2564
a_id.f54      0.23      0.06     0.11     0.35 1.00    10638     2764
a_id.f55      1.74      0.06     1.63     1.86 1.00    10136     2617
a_id.f56     -0.21      0.09    -0.37    -0.04 1.00     9053     2606
a_id.f58     -1.44      0.08    -1.60    -1.28 1.00     9842     2312
a_id.f60      0.34      0.06     0.22     0.45 1.00    10447     2865
a_id.f64      0.88      0.07     0.75     1.01 1.00    10141     2394
a_id.f65      0.78      0.07     0.64     0.91 1.00    11076     2488
a_id.f67     -0.08      0.06    -0.21     0.04 1.00    11294     2590
a_id.f68      0.04      0.06    -0.09     0.16 1.00     9932     2665
a_id.f69     -0.81      0.06    -0.93    -0.69 1.00     8488     2789
a_id.f70      0.21      0.06     0.08     0.33 1.00     8768     2299
a_id.f72     -1.21      0.06    -1.32    -1.10 1.00     9975     2808
a_id.f73     -0.12      0.06    -0.24    -0.00 1.00    10602     2661
a_id.f75      0.01      0.06    -0.11     0.12 1.00     9736     2519
a_id.f76     -0.78      0.12    -1.02    -0.53 1.00     9026     2475
a_id.f80     -0.94      0.10    -1.14    -0.75 1.00    10759     2564
a_id.f81      0.29      0.06     0.17     0.41 1.00    10492     2699
a_id.f82      0.48      0.06     0.36     0.59 1.01     8959     2513
a_id.f83     -1.63      0.06    -1.74    -1.51 1.00    10706     2766
a_id.f85     -0.48      0.06    -0.61    -0.36 1.00     8873     2712
a_id.f86     -0.24      0.08    -0.41    -0.08 1.00     9621     2851
a_id.f88      1.10      0.07     0.95     1.24 1.00    11458     2660
a_id.f91     -0.23      0.06    -0.35    -0.12 1.00    12793     2750
a_id.f92     -0.39      0.06    -0.51    -0.27 1.00     7678     2507
a_id.f94      0.02      0.06    -0.10     0.13 1.00     7897     2806
a_id.f95     -0.01      0.07    -0.15     0.13 1.00    11740     2820
a_id.f96      0.34      0.06     0.22     0.46 1.00     8942     2517
a_id.f97     -0.50      0.06    -0.62    -0.38 1.00    11475     2683
a_id.f98     -0.03      0.06    -0.15     0.09 1.00    12370     2673
a_id.f99     -0.46      0.07    -0.60    -0.33 1.00     9260     2712
a_id.f100    -1.44      0.06    -1.56    -1.31 1.00    12439     2672
a_id.f101    -0.67      0.06    -0.79    -0.56 1.00    10542     2914
a_id.f102    -1.31      0.06    -1.43    -1.19 1.00    10056     2485
a_id.f104     1.34      0.06     1.22     1.46 1.00    11217     2417
a_id.f105    -1.10      0.06    -1.22    -0.97 1.00    13252     2519
a_id.f106     0.49      0.06     0.38     0.61 1.00     7821     2608
a_id.f107     0.57      0.06     0.45     0.69 1.00     9225     2433
a_id.f109     1.66      0.06     1.55     1.78 1.00     9725     2928
a_id.f110    -0.05      0.06    -0.17     0.07 1.00     8123     2753
a_id.f111     0.55      0.06     0.44     0.67 1.00     7923     2988
a_id.f112    -0.58      0.15    -0.87    -0.29 1.00     9583     2744
a_id.f113     0.24      0.06     0.12     0.35 1.00    10205     2840
a_id.f114    -0.70      0.07    -0.84    -0.56 1.00     9191     2473
a_id.f115     0.27      0.08     0.10     0.43 1.01    10034     2335
a_id.f116     0.04      0.06    -0.08     0.16 1.00     8985     2855
a_id.f117     0.45      0.12     0.20     0.69 1.00    10768     2864
a_id.f119     0.07      0.12    -0.16     0.29 1.00     8497     2771
a_id.f120    -0.21      0.08    -0.36    -0.05 1.00     9532     2146
a_id.f121     1.22      0.06     1.10     1.34 1.00    12291     2996
a_id.f123    -0.18      0.06    -0.30    -0.06 1.00    10283     2723
a_id.f125    -0.26      0.07    -0.39    -0.13 1.00     9338     2875
a_id.f127    -1.00      0.06    -1.12    -0.88 1.00    11082     2563
a_id.f128    -1.11      0.06    -1.23    -0.98 1.00    10116     2529
a_id.f129     0.57      0.10     0.39     0.75 1.00     8745     3141
a_id.f132     0.10      0.06    -0.02     0.22 1.00    10010     2707
a_id.f133     0.02      0.06    -0.11     0.14 1.00     8711     2583
a_id.f134     0.21      0.07     0.08     0.34 1.00    10237     2613
a_id.f135     1.03      0.21     0.62     1.44 1.00    11717     2849
a_id.f136    -0.18      0.07    -0.31    -0.05 1.00    11568     2788
a_id.f137     1.66      0.08     1.51     1.83 1.01    10422     2345
a_id.f138    -1.68      0.07    -1.81    -1.55 1.00    10819     2759
a_id.f139     0.14      0.07    -0.00     0.27 1.00     9378     2591
a_id.f140     0.06      0.07    -0.07     0.19 1.00    11997     2515
a_id.f142     1.34      0.07     1.20     1.48 1.00     9723     2380
a_id.f143    -1.15      0.07    -1.28    -1.02 1.00     9321     2730
a_id.f144    -0.26      0.12    -0.50    -0.02 1.00     9379     2805
a_id.f145    -1.34      0.08    -1.50    -1.17 1.00    10669     2879
a_id.f147     0.01      0.09    -0.17     0.20 1.00    10415     2557
a_id.f148    -0.05      0.07    -0.19     0.08 1.00     9780     2429
a_id.f149     0.90      0.07     0.76     1.04 1.00    12587     2651
a_id.f150     1.55      0.07     1.42     1.69 1.00    13210     2557
a_id.f151     0.69      0.07     0.55     0.83 1.00    13447     2697
a_id.f152    -0.83      0.07    -0.97    -0.68 1.00     9650     3028
a_id.f153    -1.35      0.07    -1.50    -1.21 1.00     9672     2439
a_id.f154    -0.77      0.08    -0.92    -0.63 1.00    10363     2583
a_id.f155    -0.39      0.08    -0.54    -0.24 1.00     9490     2736
a_id.f157     0.76      0.07     0.62     0.91 1.00    12031     2628
a_id.f158     0.76      0.08     0.61     0.92 1.00    11619     2615
a_id.f160    -0.17      0.08    -0.33    -0.00 1.00     8099     2270
a_id.f162     0.36      0.08     0.20     0.52 1.00    10579     2518
a_id.f163    -0.24      0.07    -0.38    -0.09 1.00     7979     2688
a_id.f164    -1.39      0.07    -1.53    -1.25 1.00    10235     2906
a_id.f165     1.12      0.07     0.98     1.26 1.00    11826     2540
a_id.f166     0.04      0.07    -0.10     0.19 1.00     9628     2716
a_id.f167    -0.46      0.07    -0.61    -0.32 1.00    11116     2846
a_id.f168    -0.26      0.07    -0.40    -0.12 1.00     8970     2591
a_id.f169     0.36      0.08     0.21     0.52 1.01     9863     2360
a_id.f170    -1.38      0.08    -1.53    -1.23 1.00    10634     2948
a_id.f171    -0.52      0.07    -0.67    -0.37 1.00    10290     2607
a_id.f173    -0.01      0.09    -0.18     0.16 1.00     8169     2938
a_id.f174    -0.12      0.09    -0.29     0.04 1.00    12217     2605
a_id.f176     0.43      0.10     0.22     0.63 1.00     8486     2533
a_id.f177    -0.04      0.10    -0.23     0.15 1.01     8429     2566
a_id.f178    -0.19      0.08    -0.34    -0.04 1.00    12379     2182
a_id.f181     0.68      0.08     0.52     0.84 1.00    10100     2429
a_id.f182     0.40      0.08     0.24     0.56 1.00    10363     2484
a_id.f183     0.61      0.08     0.46     0.76 1.00    11405     2813
a_id.f184    -1.61      0.08    -1.77    -1.45 1.00    10245     2776
a_id.f185     0.25      0.08     0.10     0.40 1.00     9114     2973
a_id.f186    -0.62      0.09    -0.80    -0.45 1.00     8718     2540
a_id.f188    -0.20      0.08    -0.36    -0.04 1.00    10469     2604
a_id.f189     0.26      0.10     0.07     0.45 1.00    11539     2702
a_id.f190    -1.72      0.08    -1.89    -1.56 1.00     8338     2445
a_id.f191    -0.06      0.09    -0.22     0.11 1.00    11183     2658
a_id.f194     0.02      0.09    -0.16     0.20 1.00     9761     2193
a_id.f195     1.43      0.08     1.26     1.60 1.00     8897     3002
a_id.f196     0.02      0.10    -0.17     0.20 1.00     7595     2485
a_id.f197     0.64      0.09     0.47     0.81 1.00    11572     2250
a_id.f198    -1.65      0.09    -1.82    -1.46 1.00     8803     3035
a_id.f199     0.46      0.09     0.28     0.63 1.00     9935     2806
a_id.f201     0.68      0.08     0.53     0.82 1.00     8205     2792
a_id.f202    -0.47      0.09    -0.64    -0.29 1.00    10762     3021
a_id.f203    -1.10      0.10    -1.30    -0.90 1.00    10729     2514
a_id.f204    -1.20      0.10    -1.40    -0.99 1.00     9415     2508
a_id.f205    -0.00      0.10    -0.19     0.19 1.00     9123     2452
a_id.f206    -0.18      0.09    -0.35    -0.01 1.00     9946     2997
a_id.f207    -0.41      0.10    -0.60    -0.22 1.00    12955     2891
a_id.f208     0.35      0.10     0.15     0.55 1.00    12311     2147
a_id.f209    -0.47      0.10    -0.67    -0.28 1.00     9681     2707
a_id.f211    -1.27      0.08    -1.43    -1.11 1.00     9841     2794
a_id.f212     0.12      0.10    -0.06     0.32 1.00    10195     2916
a_id.f213     0.17      0.12    -0.07     0.41 1.00    11124     2958
a_id.f214     0.42      0.11     0.21     0.64 1.00     8942     2451
a_id.f216    -0.06      0.13    -0.32     0.19 1.00    11152     2525
a_id.f217     0.03      0.10    -0.17     0.23 1.00    10467     2976
a_id.f219     0.07      0.10    -0.14     0.27 1.00     9617     2229
a_id.f220     1.01      0.11     0.79     1.22 1.00     9917     2500
a_id.f221    -0.26      0.11    -0.48    -0.04 1.00    11306     2733
a_id.f222     0.16      0.10    -0.04     0.37 1.00     8358     2780
a_id.f223    -0.12      0.12    -0.36     0.12 1.00     9828     2595
a_id.f224    -0.24      0.12    -0.48    -0.01 1.00     9876     2741
a_id.f225    -0.09      0.12    -0.32     0.14 1.00     9183     2742
a_id.f226    -0.74      0.14    -1.02    -0.46 1.00    11079     2789
a_id.f227    -0.06      0.15    -0.37     0.25 1.00     8192     2470
a_id.f228     0.44      0.14     0.17     0.72 1.00     9356     2342
a_id.f229     0.52      0.14     0.24     0.80 1.00    10801     2699
a_id.f230    -0.42      0.14    -0.71    -0.14 1.00    10968     2590
a_id.f231     0.59      0.16     0.28     0.90 1.00    10388     2354
a_id.f232     0.60      0.12     0.37     0.83 1.00    10335     2461
a_id.f233     0.39      0.12     0.14     0.63 1.01     9310     2710
a_id.f234     0.90      0.17     0.55     1.23 1.00     9510     2748
a_id.f236    -1.32      0.14    -1.58    -1.06 1.00     9147     2398
a_id.f237    -0.63      0.15    -0.91    -0.34 1.00     9507     2838
a_id.f238     1.04      0.15     0.76     1.32 1.00    10593     2814
a_id.f239     0.55      0.18     0.20     0.90 1.01    10202     2706

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.60      0.00     0.59     0.61 1.00     7117     2734

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

This is inefficient, in that the model treat each person as entirely separate. Let’s try a partial pooling model.

\[\begin{align*} \text{PA.std}_i &\sim \text{Normal}(\mu_i, \sigma) \\ \mu_i &= \alpha_{\text{id}[i]} \\ \alpha_j &\sim \text{Normal}(\bar{\alpha}, \sigma_{\alpha}) \text{ for j in 1...239}\\ \bar{\alpha} &\sim \text{Normal}(0, 1.5)\\ \sigma_{\alpha} &\sim \text{Exponential}(1) \\ \sigma \sim &\text{Exponential}(1) \\ \end{align*}\]

m2 <- brm(
  data=d,
  family=gaussian,
  PA.std ~ 1 + (1 | id), 
  prior = c( prior(normal(0, 1.5), class=Intercept),
             prior(exponential(1), class=sd),
             prior(exponential(1), class=sigma)),
  iter=4000, warmup=1000, chains=4, cores=4, seed=9, #running this longer so it mixes
  file=here("files/models/m71.2")
)
m2
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: PA.std ~ 1 + (1 | id) 
   Data: d (Number of observations: 13033) 
  Draws: 4 chains, each with iter = 4000; warmup = 1000; thin = 1;
         total post-warmup draws = 12000

Multilevel Hyperparameters:
~id (Number of levels: 193) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.80      0.04     0.73     0.89 1.01      411      941

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept    -0.01      0.06    -0.12     0.10 1.01      137      292

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.60      0.00     0.59     0.61 1.00    15176     9021

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

(BTW, I have some poor mixing here)

mcmc_rank_overlay(m2, pars = "sd_id__Intercept")

Options:

  • run with more iterations per chain.
  • increase adapt_delta to .99.

(Probably do both).

how many parameters does each model have?

get_variables(m1) 
  [1] "b_a_id.f1"     "b_a_id.f2"     "b_a_id.f3"     "b_a_id.f4"    
  [5] "b_a_id.f5"     "b_a_id.f6"     "b_a_id.f7"     "b_a_id.f8"    
  [9] "b_a_id.f9"     "b_a_id.f10"    "b_a_id.f11"    "b_a_id.f12"   
 [13] "b_a_id.f13"    "b_a_id.f14"    "b_a_id.f15"    "b_a_id.f16"   
 [17] "b_a_id.f17"    "b_a_id.f20"    "b_a_id.f23"    "b_a_id.f24"   
 [21] "b_a_id.f25"    "b_a_id.f26"    "b_a_id.f27"    "b_a_id.f28"   
 [25] "b_a_id.f29"    "b_a_id.f30"    "b_a_id.f31"    "b_a_id.f32"   
 [29] "b_a_id.f33"    "b_a_id.f34"    "b_a_id.f35"    "b_a_id.f36"   
 [33] "b_a_id.f37"    "b_a_id.f38"    "b_a_id.f39"    "b_a_id.f40"   
 [37] "b_a_id.f41"    "b_a_id.f42"    "b_a_id.f43"    "b_a_id.f44"   
 [41] "b_a_id.f46"    "b_a_id.f47"    "b_a_id.f48"    "b_a_id.f49"   
 [45] "b_a_id.f50"    "b_a_id.f51"    "b_a_id.f52"    "b_a_id.f53"   
 [49] "b_a_id.f54"    "b_a_id.f55"    "b_a_id.f56"    "b_a_id.f58"   
 [53] "b_a_id.f60"    "b_a_id.f64"    "b_a_id.f65"    "b_a_id.f67"   
 [57] "b_a_id.f68"    "b_a_id.f69"    "b_a_id.f70"    "b_a_id.f72"   
 [61] "b_a_id.f73"    "b_a_id.f75"    "b_a_id.f76"    "b_a_id.f80"   
 [65] "b_a_id.f81"    "b_a_id.f82"    "b_a_id.f83"    "b_a_id.f85"   
 [69] "b_a_id.f86"    "b_a_id.f88"    "b_a_id.f91"    "b_a_id.f92"   
 [73] "b_a_id.f94"    "b_a_id.f95"    "b_a_id.f96"    "b_a_id.f97"   
 [77] "b_a_id.f98"    "b_a_id.f99"    "b_a_id.f100"   "b_a_id.f101"  
 [81] "b_a_id.f102"   "b_a_id.f104"   "b_a_id.f105"   "b_a_id.f106"  
 [85] "b_a_id.f107"   "b_a_id.f109"   "b_a_id.f110"   "b_a_id.f111"  
 [89] "b_a_id.f112"   "b_a_id.f113"   "b_a_id.f114"   "b_a_id.f115"  
 [93] "b_a_id.f116"   "b_a_id.f117"   "b_a_id.f119"   "b_a_id.f120"  
 [97] "b_a_id.f121"   "b_a_id.f123"   "b_a_id.f125"   "b_a_id.f127"  
[101] "b_a_id.f128"   "b_a_id.f129"   "b_a_id.f132"   "b_a_id.f133"  
[105] "b_a_id.f134"   "b_a_id.f135"   "b_a_id.f136"   "b_a_id.f137"  
[109] "b_a_id.f138"   "b_a_id.f139"   "b_a_id.f140"   "b_a_id.f142"  
[113] "b_a_id.f143"   "b_a_id.f144"   "b_a_id.f145"   "b_a_id.f147"  
[117] "b_a_id.f148"   "b_a_id.f149"   "b_a_id.f150"   "b_a_id.f151"  
[121] "b_a_id.f152"   "b_a_id.f153"   "b_a_id.f154"   "b_a_id.f155"  
[125] "b_a_id.f157"   "b_a_id.f158"   "b_a_id.f160"   "b_a_id.f162"  
[129] "b_a_id.f163"   "b_a_id.f164"   "b_a_id.f165"   "b_a_id.f166"  
[133] "b_a_id.f167"   "b_a_id.f168"   "b_a_id.f169"   "b_a_id.f170"  
[137] "b_a_id.f171"   "b_a_id.f173"   "b_a_id.f174"   "b_a_id.f176"  
[141] "b_a_id.f177"   "b_a_id.f178"   "b_a_id.f181"   "b_a_id.f182"  
[145] "b_a_id.f183"   "b_a_id.f184"   "b_a_id.f185"   "b_a_id.f186"  
[149] "b_a_id.f188"   "b_a_id.f189"   "b_a_id.f190"   "b_a_id.f191"  
[153] "b_a_id.f194"   "b_a_id.f195"   "b_a_id.f196"   "b_a_id.f197"  
[157] "b_a_id.f198"   "b_a_id.f199"   "b_a_id.f201"   "b_a_id.f202"  
[161] "b_a_id.f203"   "b_a_id.f204"   "b_a_id.f205"   "b_a_id.f206"  
[165] "b_a_id.f207"   "b_a_id.f208"   "b_a_id.f209"   "b_a_id.f211"  
[169] "b_a_id.f212"   "b_a_id.f213"   "b_a_id.f214"   "b_a_id.f216"  
[173] "b_a_id.f217"   "b_a_id.f219"   "b_a_id.f220"   "b_a_id.f221"  
[177] "b_a_id.f222"   "b_a_id.f223"   "b_a_id.f224"   "b_a_id.f225"  
[181] "b_a_id.f226"   "b_a_id.f227"   "b_a_id.f228"   "b_a_id.f229"  
[185] "b_a_id.f230"   "b_a_id.f231"   "b_a_id.f232"   "b_a_id.f233"  
[189] "b_a_id.f234"   "b_a_id.f236"   "b_a_id.f237"   "b_a_id.f238"  
[193] "b_a_id.f239"   "sigma"         "lprior"        "lp__"         
[197] "accept_stat__" "stepsize__"    "treedepth__"   "n_leapfrog__" 
[201] "divergent__"   "energy__"

how many parameters does each model have?

get_variables(m1) %>% length()
[1] 202
get_variables(m2) %>% length()
[1] 205

What additional parameters?

m1 has a unique intercept for each participant and a standard deviation of scores (1 \(\sigma\)).

m2 is estimating all of that plus a grand mean intercept and the variability of means (\(\sigma_M\)).

(what’s the extra one? brms lists the intercept twice. *shrug emoji*)

And yet!

m1 <- add_criterion(m1, criterion = "loo")
m2 <- add_criterion(m2, criterion = "loo")

loo_compare(m1, m2) %>% print(simplify=F)
   elpd_diff se_diff  elpd_loo se_elpd_loo p_loo    se_p_loo looic    se_looic
m2      0.0       0.0 -11955.9    110.5       201.5      4.4  23911.8    221.0
m1     -0.1       0.9 -11956.0    110.5       202.1      4.4  23911.9    221.1

Let’s visualize the differences in these.

Code 
nd1 = distinct(d, id.f)
post1 = epred_draws(m1, nd1)
nd2 = distinct(d, id)
post2 = epred_draws(m2, nd2)
p1 = post1 %>% 
  ggplot( aes(y=.epred, x=id.f) ) +
  stat_gradientinterval() +
  scale_x_discrete(labels=NULL, breaks=NULL) +
  labs(x="id", y="PA.std", title = "no pooling")

p2 = post2 %>% 
  mutate(id=as.factor(id)) %>% 
  ggplot( aes(y=.epred, x=id) ) +
  stat_gradientinterval() +
  scale_x_discrete(labels=NULL, breaks=NULL) +
  labs(x="id", y="PA.std", title = "partial pooling")

p1 / p2
Code 
means1 = post1 %>% 
  mean_qi(.epred)
means2 = post2 %>% 
  mean_qi(.epred) %>% 
  mutate(id=as.factor(id))

means1 %>% 
  ggplot( aes(x=id.f, y=.epred)) +
  geom_hline( aes(yintercept=mean(.epred)),
              linetype="dashed") +
  geom_point( aes(color="no pooling") ) +
  geom_point( aes(x=id, color="partial pooling"),
              data=means2,
              size=2,
              alpha=.4) +
  scale_color_manual( values=c("#e07a5f", "#1c5253") ) +
  scale_x_discrete(breaks=NULL) +
  labs(x="id", y="PA.std")+
  theme(legend.position = "top")

The amount of pooling is affected by the number of data points – lots of data for each cluster = less pooling.

Code 
d_small = d %>% 
  with_groups(id, 
              sample_n,
              size = 3)
m1s <- brm(
  data=d_small,
  family=gaussian,
  bf(PA.std ~ 0 + a,
     a ~ 0 + id.f,
     nl = TRUE),
  prior = c( prior(normal(0, 1.5), class=b, nlpar=a),
             prior(exponential(1), class=sigma)),
  iter=2000, warmup=1000, chains=4, cores=4, seed=9,
  file=here("files/models/m71.1s")
)

m2s <- brm(
  data=d_small,
  family=gaussian,
  PA.std ~ 1 + (1 | id), 
  prior = c( prior(normal(0, 1.5), class=Intercept),
             prior(exponential(1), class=sd),
             prior(exponential(1), class=sigma)),
  iter=4000, warmup=1000, chains=4, cores=4, seed=9, #running this longer so it mixes
  file=here("files/models/m71.2s")
)

nd1 = distinct(d_small, id.f)
post1 = epred_draws(m1s, nd1)
nd2 = distinct(d_small, id)
post2 = epred_draws(m2s, nd2)

means1 = post1 %>% 
  mean_qi(.epred)
means2 = post2 %>% 
  mean_qi(.epred) %>% 
  mutate(id=as.factor(id))

means1 %>% 
  ggplot( aes(x=id.f, y=.epred)) +
  geom_hline( aes(yintercept=mean(.epred)),
              linetype="dashed") +
  geom_point( aes(color="no pooling") ) +
  geom_point( aes(x=id, color="partial pooling"),
              data=means2,
              size=2,
              alpha=.4) +
  scale_color_manual( values=c("#e07a5f", "#1c5253") ) +
  scale_x_discrete(breaks=NULL) +
  labs(x="id", y="PA.std", title="3 observations per person")+
  theme(legend.position = "bottom")

\[\begin{align*} \text{PA.std}_{ij} &\sim \text{Normal}(\mu_i, \sigma) \\ \mu_i &= \alpha_{id[i]} \\ \alpha_{id[i]} &\sim \text{Normal}(\bar{\alpha}, 1) \text{ for i in } 1...239\\ \bar{\alpha} &\sim \text{Normal}(.50,.25) \\ \sigma &\sim \text{Exponential}(1) \\ \end{align*}\]

Another way to write out the expected value of \(\mu_j\) is:

\[ \mu_j = \bar{\alpha} + \alpha_{i[ID]} \]

Where \(\bar{\alpha}\) is the mean of the means and \(\alpha_{i[ID]}\) is the deviation of each person from the grand mean. There will be one \(\bar{\alpha}\) for the entire sample, but a different \(\alpha_{i[ID]}\) for each person. Because there are multiple values for each person, we can calculate the variability of those \(\alpha\)’s, as well as the residual variability within each person.

\[\begin{align*} \text{PA.std}_{ij} &\sim \text{Normal}(\mu_j, \sigma) \\ \mu_i &= \bar{\alpha} + \alpha_{id[j]} \\ \bar{\alpha} &\sim \text{Normal}(.50,.25) \\ \sigma &\sim \text{Exponential}(1) \\ \sigma_{\alpha} &\sim \text{Exponential}(1) \\ \end{align*}\]

In this reparameterization, we no longer need a prior for each of the person-means. That’s because they must have a mean of 0. Instead, we only need to estimate the variability \((\sigma_{\alpha})\) of these deviations.

In other words, this model is analogous to an ANOVA, in which we calculate both within- and between-group (aka person) variability. Therefore, we’ll have a standard deviation for both levels.

magical tidybayes

spread_draws(m2, b_Intercept, r_id[id, term])
# A tibble: 2,316,000 × 7
# Groups:   id, term [193]
   .chain .iteration .draw b_Intercept    id term         r_id
    <int>      <int> <int>       <dbl> <int> <chr>       <dbl>
 1      1          1     1     0.00562     1 Intercept -0.0317
 2      1          1     1     0.00562     2 Intercept  1.85  
 3      1          1     1     0.00562     3 Intercept  0.0469
 4      1          1     1     0.00562     4 Intercept -0.332 
 5      1          1     1     0.00562     5 Intercept  1.24  
 6      1          1     1     0.00562     6 Intercept -0.546 
 7      1          1     1     0.00562     7 Intercept  0.202 
 8      1          1     1     0.00562     8 Intercept  2.09  
 9      1          1     1     0.00562     9 Intercept  1.01  
10      1          1     1     0.00562    10 Intercept  0.462 
# ℹ 2,315,990 more rows
Code 
spread_draws(m2, b_Intercept, r_id[id, term]) %>% 
  mutate(p_mean = b_Intercept + r_id) %>% 
  filter(id %in% sample_id) %>% 
  ggplot(aes(x = id, y = p_mean)) + 
  geom_hline( aes(yintercept = mean(d$PA.std)),
              linetype="dashed") +
  stat_halfeye(alpha =.5)

In this figure, we can visualize three approaches to pooling:

  • NO POOLING: each person’s mean is calculated from only their own data, there is no sharing of information. (black dots)
  • COMPLETE POOLING: any differences between people are just random noise. The mean of the whole group is the only estimate and is assumed to apply equally to each person (dashed line).
  • PARTIAL POOLING: information is shared across participants, but they are assumed to be different. We regularize estimates of person-specific means towards the grand mean. The more information we have on a person, the less their estimate is regularized.
Code 
actual_means = d %>% 
  with_groups(id, summarise, p_mean = mean(PA.std)) %>% 
  mutate()

spread_draws(m2, b_Intercept, r_id[id, term]) %>% 
  mutate(p_mean = b_Intercept + r_id) %>% 
  mean_qi(p_mean) %>% 
  filter(id %in% sample_id) %>% 
  ggplot( aes(x=id, y=p_mean) ) +
  geom_hline( aes(yintercept = mean(d$PA.std)),
              linetype="dashed") +
  geom_point( size=2, color="#e07a5f") + 
  geom_point( data=filter(actual_means, 
                          id %in% sample_id)) +
  labs(x="id",y="PA.std")

writing our mulitilevel model

It’s common to write out mulitlevel models using formulas for the different levels. Level 1 is the level of your outcome, or the thing that repeats. Level 2 is the level of your groups.

\[\begin{align*} \text{Level 1} &\\ \text{PA.std}_{ij} &= \alpha_j + \epsilon_{ij} \\ \text{Level 2} &\\ \alpha_j &= \gamma_0 + U_j \end{align*}\]

Some refer to \(U_j\) as a “random” or “varying” effect because it varies across groups. \(\gamma_0\) is therefore a “fixed” or “non-varying” effect because it applies to all groups.

drawing from the posterior

We’ve already seen how we can use spread_draws to get our parameter estimates.

m2 %>% spread_draws(b_Intercept, sigma, sd_id__Intercept) %>% head
# A tibble: 6 × 6
  .chain .iteration .draw b_Intercept sigma sd_id__Intercept
   <int>      <int> <int>       <dbl> <dbl>            <dbl>
1      1          1     1   0.00562   0.603            0.877
2      1          2     2  -0.0139    0.604            0.873
3      1          3     3  -0.00327   0.599            0.830
4      1          4     4  -0.0000504 0.597            0.819
5      1          5     5  -0.0114    0.599            0.818
6      1          6     6  -0.0238    0.599            0.818
m2 %>% spread_draws(r_id[id, term]) %>% head
# A tibble: 6 × 6
# Groups:   id, term [1]
     id term         r_id .chain .iteration .draw
  <int> <chr>       <dbl>  <int>      <int> <int>
1     1 Intercept -0.0317      1          1     1
2     1 Intercept -0.0562      1          2     2
3     1 Intercept  0.228       1          3     3
4     1 Intercept  0.218       1          4     4
5     1 Intercept  0.0904      1          5     5
6     1 Intercept  0.159       1          6     6

We can get expected values (means) for individuals in our sample.

set.seed(9)
sample_id = sample(unique(d$id), replace=F, size=3)

nd = data.frame(id=sample_id, id.f = sample_id)

add_epred_draws(newdata=nd, object=m2) %>% 
  mean_qi()
# A tibble: 3 × 9
     id  id.f  .row .epred .lower .upper .width .point .interval
  <int> <int> <int>  <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1     6     6     3 -0.525 -0.664 -0.386   0.95 mean   qi       
2    60    60     2  0.335  0.213  0.456   0.95 mean   qi       
3   232   232     1  0.591  0.358  0.831   0.95 mean   qi

We can get predicted values (scores) for individuals in our sample.

add_predicted_draws(newdata=nd, object=m2) %>% 
  mean_qi()
# A tibble: 3 × 9
     id  id.f  .row .prediction .lower .upper .width .point .interval
  <int> <int> <int>       <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1     6     6     3      -0.531 -1.72   0.644   0.95 mean   qi       
2    60    60     2       0.325 -0.885  1.52    0.95 mean   qi       
3   232   232     1       0.583 -0.614  1.76    0.95 mean   qi

expected values vs predicted values

Code 
expected = add_epred_draws(newdata=nd, object=m2)
predicted = add_predicted_draws(newdata=nd, object=m2)

full_join(expected, predicted) %>% 
  filter(.draw <= 100) %>% 
  ungroup() %>% 
  ggplot( aes(x=id, y=.epred)) +
  stat_halfeye( aes(fill="expected"),
                alpha=.3) +
  stat_halfeye( aes(y =.prediction, fill="predicted"), 
                alpha=.3 ) +
  scale_fill_manual( values=c("#e07a5f","#1c5253") ) +
  labs(x="id", y="PA.std", fill="draw") +
  scale_x_continuous(breaks=sample_id) +
  theme(legend.position = "top")

Finally, we can get expected values for new individuals

Code 
new_id = max(d$id) + seq(1:3)

nd = data.frame(id=new_id)

add_epred_draws(newdata=nd, object=m2, 
                allow_new_levels=T) %>% 
  ggplot( aes(x=id, y=.epred) ) +
  stat_halfeye()

adding covariates

Let’s add time to our model. Because time varies (can change) within person from assessment to assessment, this is a Level 1 variable. Note that I have NOT added a varying component to Level 2 – in other words, I’m stating that the effect of time on positive affect is fixed or identical across participants.

\[\begin{align*} \text{Level 1} &\\ \text{PA.std}_{ij} &= \alpha_i + \beta_i(\text{day}_{ij}) + \epsilon_{ij} \\ \text{Level 2} &\\ \alpha_j &= \gamma_0 + U_{0j} \\ \beta_j &= \gamma_1 \\ \end{align*}\]

m3 <- brm(
  data=d,
  family=gaussian,
  PA.std ~ 1 + day + (1 | id), 
  prior = c( prior(normal(.50, .25), class=Intercept),
             prior(normal(0, 1), class=b), #new prior for new term
             prior(exponential(1), class=sd),
             prior(exponential(1), class=sigma)),
  iter=2000, warmup=1000, chains=4, cores=4, seed=9,
  file=here("files/models/m71.3")
)
m3
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: PA.std ~ 1 + day + (1 | id) 
   Data: d (Number of observations: 13033) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Multilevel Hyperparameters:
~id (Number of levels: 193) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.80      0.04     0.73     0.88 1.04      110      285

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     0.13      0.05     0.02     0.23 1.08       47      188
day          -0.00      0.00    -0.00    -0.00 1.00     4550     3272

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.60      0.00     0.59     0.60 1.00     6598     2554

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

There are different intercepts for each participant, but not different slopes.

get_variables(m3)
  [1] "b_Intercept"         "b_day"               "sd_id__Intercept"   
  [4] "sigma"               "Intercept"           "r_id[1,Intercept]"  
  [7] "r_id[2,Intercept]"   "r_id[3,Intercept]"   "r_id[4,Intercept]"  
 [10] "r_id[5,Intercept]"   "r_id[6,Intercept]"   "r_id[7,Intercept]"  
 [13] "r_id[8,Intercept]"   "r_id[9,Intercept]"   "r_id[10,Intercept]" 
 [16] "r_id[11,Intercept]"  "r_id[12,Intercept]"  "r_id[13,Intercept]" 
 [19] "r_id[14,Intercept]"  "r_id[15,Intercept]"  "r_id[16,Intercept]" 
 [22] "r_id[17,Intercept]"  "r_id[20,Intercept]"  "r_id[23,Intercept]" 
 [25] "r_id[24,Intercept]"  "r_id[25,Intercept]"  "r_id[26,Intercept]" 
 [28] "r_id[27,Intercept]"  "r_id[28,Intercept]"  "r_id[29,Intercept]" 
 [31] "r_id[30,Intercept]"  "r_id[31,Intercept]"  "r_id[32,Intercept]" 
 [34] "r_id[33,Intercept]"  "r_id[34,Intercept]"  "r_id[35,Intercept]" 
 [37] "r_id[36,Intercept]"  "r_id[37,Intercept]"  "r_id[38,Intercept]" 
 [40] "r_id[39,Intercept]"  "r_id[40,Intercept]"  "r_id[41,Intercept]" 
 [43] "r_id[42,Intercept]"  "r_id[43,Intercept]"  "r_id[44,Intercept]" 
 [46] "r_id[46,Intercept]"  "r_id[47,Intercept]"  "r_id[48,Intercept]" 
 [49] "r_id[49,Intercept]"  "r_id[50,Intercept]"  "r_id[51,Intercept]" 
 [52] "r_id[52,Intercept]"  "r_id[53,Intercept]"  "r_id[54,Intercept]" 
 [55] "r_id[55,Intercept]"  "r_id[56,Intercept]"  "r_id[58,Intercept]" 
 [58] "r_id[60,Intercept]"  "r_id[64,Intercept]"  "r_id[65,Intercept]" 
 [61] "r_id[67,Intercept]"  "r_id[68,Intercept]"  "r_id[69,Intercept]" 
 [64] "r_id[70,Intercept]"  "r_id[72,Intercept]"  "r_id[73,Intercept]" 
 [67] "r_id[75,Intercept]"  "r_id[76,Intercept]"  "r_id[80,Intercept]" 
 [70] "r_id[81,Intercept]"  "r_id[82,Intercept]"  "r_id[83,Intercept]" 
 [73] "r_id[85,Intercept]"  "r_id[86,Intercept]"  "r_id[88,Intercept]" 
 [76] "r_id[91,Intercept]"  "r_id[92,Intercept]"  "r_id[94,Intercept]" 
 [79] "r_id[95,Intercept]"  "r_id[96,Intercept]"  "r_id[97,Intercept]" 
 [82] "r_id[98,Intercept]"  "r_id[99,Intercept]"  "r_id[100,Intercept]"
 [85] "r_id[101,Intercept]" "r_id[102,Intercept]" "r_id[104,Intercept]"
 [88] "r_id[105,Intercept]" "r_id[106,Intercept]" "r_id[107,Intercept]"
 [91] "r_id[109,Intercept]" "r_id[110,Intercept]" "r_id[111,Intercept]"
 [94] "r_id[112,Intercept]" "r_id[113,Intercept]" "r_id[114,Intercept]"
 [97] "r_id[115,Intercept]" "r_id[116,Intercept]" "r_id[117,Intercept]"
[100] "r_id[119,Intercept]" "r_id[120,Intercept]" "r_id[121,Intercept]"
[103] "r_id[123,Intercept]" "r_id[125,Intercept]" "r_id[127,Intercept]"
[106] "r_id[128,Intercept]" "r_id[129,Intercept]" "r_id[132,Intercept]"
[109] "r_id[133,Intercept]" "r_id[134,Intercept]" "r_id[135,Intercept]"
[112] "r_id[136,Intercept]" "r_id[137,Intercept]" "r_id[138,Intercept]"
[115] "r_id[139,Intercept]" "r_id[140,Intercept]" "r_id[142,Intercept]"
[118] "r_id[143,Intercept]" "r_id[144,Intercept]" "r_id[145,Intercept]"
[121] "r_id[147,Intercept]" "r_id[148,Intercept]" "r_id[149,Intercept]"
[124] "r_id[150,Intercept]" "r_id[151,Intercept]" "r_id[152,Intercept]"
[127] "r_id[153,Intercept]" "r_id[154,Intercept]" "r_id[155,Intercept]"
[130] "r_id[157,Intercept]" "r_id[158,Intercept]" "r_id[160,Intercept]"
[133] "r_id[162,Intercept]" "r_id[163,Intercept]" "r_id[164,Intercept]"
[136] "r_id[165,Intercept]" "r_id[166,Intercept]" "r_id[167,Intercept]"
[139] "r_id[168,Intercept]" "r_id[169,Intercept]" "r_id[170,Intercept]"
[142] "r_id[171,Intercept]" "r_id[173,Intercept]" "r_id[174,Intercept]"
[145] "r_id[176,Intercept]" "r_id[177,Intercept]" "r_id[178,Intercept]"
[148] "r_id[181,Intercept]" "r_id[182,Intercept]" "r_id[183,Intercept]"
[151] "r_id[184,Intercept]" "r_id[185,Intercept]" "r_id[186,Intercept]"
[154] "r_id[188,Intercept]" "r_id[189,Intercept]" "r_id[190,Intercept]"
[157] "r_id[191,Intercept]" "r_id[194,Intercept]" "r_id[195,Intercept]"
[160] "r_id[196,Intercept]" "r_id[197,Intercept]" "r_id[198,Intercept]"
[163] "r_id[199,Intercept]" "r_id[201,Intercept]" "r_id[202,Intercept]"
[166] "r_id[203,Intercept]" "r_id[204,Intercept]" "r_id[205,Intercept]"
[169] "r_id[206,Intercept]" "r_id[207,Intercept]" "r_id[208,Intercept]"
[172] "r_id[209,Intercept]" "r_id[211,Intercept]" "r_id[212,Intercept]"
[175] "r_id[213,Intercept]" "r_id[214,Intercept]" "r_id[216,Intercept]"
[178] "r_id[217,Intercept]" "r_id[219,Intercept]" "r_id[220,Intercept]"
[181] "r_id[221,Intercept]" "r_id[222,Intercept]" "r_id[223,Intercept]"
[184] "r_id[224,Intercept]" "r_id[225,Intercept]" "r_id[226,Intercept]"
[187] "r_id[227,Intercept]" "r_id[228,Intercept]" "r_id[229,Intercept]"
[190] "r_id[230,Intercept]" "r_id[231,Intercept]" "r_id[232,Intercept]"
[193] "r_id[233,Intercept]" "r_id[234,Intercept]" "r_id[236,Intercept]"
[196] "r_id[237,Intercept]" "r_id[238,Intercept]" "r_id[239,Intercept]"
[199] "lprior"              "lp__"                "accept_stat__"      
[202] "stepsize__"          "treedepth__"         "n_leapfrog__"       
[205] "divergent__"         "energy__"
Code 
set.seed(9)
sample_id = sample(unique(d$id), replace=F, size = 20)
distinct(d, id, day) %>% 
  filter(id %in% sample_id) %>% 
  add_epred_draws(m3) %>% 
  mean_qi() %>% 
  ggplot( aes(x=day, y=.epred, color=as.factor(id))) +
  geom_line() +
  guides(color="none") +
  labs(x="day",y="PA.std")

level 2 covariates

We can also add covariates at Level 2, or the person-level in this case. We have a binary variable called female that refers to the gender of each person1.

\[\begin{align*} \text{Level 1} &\\ \text{PA.std}_{ij} &= \alpha_i + \beta_i(\text{day}_{ij}) + \epsilon_{ij} \\ \text{Level 2} &\\ \alpha_j &= \gamma_{0\text{female}[j]} + U_{0j} \\ \beta_j &= \gamma_{1} \\ \end{align*}\]

m4 <- brm(
  data=d,
  family=gaussian,
    PA.std ~ 0 + day + female + (1 | id), 
  prior = c( prior(normal(0, 1), class=b), 
             prior(exponential(1), class=sd),
             prior(exponential(1), class=sigma)),
  iter=5000, warmup=1000, chains=4, cores=4, seed=9,
  file=here("files/models/m71.4")
)
m4
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: PA.std ~ 0 + day + female + (1 | id) 
   Data: d (Number of observations: 13033) 
  Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 1;
         total post-warmup draws = 16000

Multilevel Hyperparameters:
~id (Number of levels: 193) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.77      0.04     0.70     0.85 1.00      632     1163

Regression Coefficients:
       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
day       -0.00      0.00    -0.00    -0.00 1.00    13839    13412
female     0.35      0.07     0.20     0.49 1.01      371      731

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.60      0.00     0.59     0.60 1.00    36313    10598

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
get_variables(m4)
  [1] "b_day"               "b_female"            "sd_id__Intercept"   
  [4] "sigma"               "r_id[1,Intercept]"   "r_id[2,Intercept]"  
  [7] "r_id[3,Intercept]"   "r_id[4,Intercept]"   "r_id[5,Intercept]"  
 [10] "r_id[6,Intercept]"   "r_id[7,Intercept]"   "r_id[8,Intercept]"  
 [13] "r_id[9,Intercept]"   "r_id[10,Intercept]"  "r_id[11,Intercept]" 
 [16] "r_id[12,Intercept]"  "r_id[13,Intercept]"  "r_id[14,Intercept]" 
 [19] "r_id[15,Intercept]"  "r_id[16,Intercept]"  "r_id[17,Intercept]" 
 [22] "r_id[20,Intercept]"  "r_id[23,Intercept]"  "r_id[24,Intercept]" 
 [25] "r_id[25,Intercept]"  "r_id[26,Intercept]"  "r_id[27,Intercept]" 
 [28] "r_id[28,Intercept]"  "r_id[29,Intercept]"  "r_id[30,Intercept]" 
 [31] "r_id[31,Intercept]"  "r_id[32,Intercept]"  "r_id[33,Intercept]" 
 [34] "r_id[34,Intercept]"  "r_id[35,Intercept]"  "r_id[36,Intercept]" 
 [37] "r_id[37,Intercept]"  "r_id[38,Intercept]"  "r_id[39,Intercept]" 
 [40] "r_id[40,Intercept]"  "r_id[41,Intercept]"  "r_id[42,Intercept]" 
 [43] "r_id[43,Intercept]"  "r_id[44,Intercept]"  "r_id[46,Intercept]" 
 [46] "r_id[47,Intercept]"  "r_id[48,Intercept]"  "r_id[49,Intercept]" 
 [49] "r_id[50,Intercept]"  "r_id[51,Intercept]"  "r_id[52,Intercept]" 
 [52] "r_id[53,Intercept]"  "r_id[54,Intercept]"  "r_id[55,Intercept]" 
 [55] "r_id[56,Intercept]"  "r_id[58,Intercept]"  "r_id[60,Intercept]" 
 [58] "r_id[64,Intercept]"  "r_id[65,Intercept]"  "r_id[67,Intercept]" 
 [61] "r_id[68,Intercept]"  "r_id[69,Intercept]"  "r_id[70,Intercept]" 
 [64] "r_id[72,Intercept]"  "r_id[73,Intercept]"  "r_id[75,Intercept]" 
 [67] "r_id[76,Intercept]"  "r_id[80,Intercept]"  "r_id[81,Intercept]" 
 [70] "r_id[82,Intercept]"  "r_id[83,Intercept]"  "r_id[85,Intercept]" 
 [73] "r_id[86,Intercept]"  "r_id[88,Intercept]"  "r_id[91,Intercept]" 
 [76] "r_id[92,Intercept]"  "r_id[94,Intercept]"  "r_id[95,Intercept]" 
 [79] "r_id[96,Intercept]"  "r_id[97,Intercept]"  "r_id[98,Intercept]" 
 [82] "r_id[99,Intercept]"  "r_id[100,Intercept]" "r_id[101,Intercept]"
 [85] "r_id[102,Intercept]" "r_id[104,Intercept]" "r_id[105,Intercept]"
 [88] "r_id[106,Intercept]" "r_id[107,Intercept]" "r_id[109,Intercept]"
 [91] "r_id[110,Intercept]" "r_id[111,Intercept]" "r_id[112,Intercept]"
 [94] "r_id[113,Intercept]" "r_id[114,Intercept]" "r_id[115,Intercept]"
 [97] "r_id[116,Intercept]" "r_id[117,Intercept]" "r_id[119,Intercept]"
[100] "r_id[120,Intercept]" "r_id[121,Intercept]" "r_id[123,Intercept]"
[103] "r_id[125,Intercept]" "r_id[127,Intercept]" "r_id[128,Intercept]"
[106] "r_id[129,Intercept]" "r_id[132,Intercept]" "r_id[133,Intercept]"
[109] "r_id[134,Intercept]" "r_id[135,Intercept]" "r_id[136,Intercept]"
[112] "r_id[137,Intercept]" "r_id[138,Intercept]" "r_id[139,Intercept]"
[115] "r_id[140,Intercept]" "r_id[142,Intercept]" "r_id[143,Intercept]"
[118] "r_id[144,Intercept]" "r_id[145,Intercept]" "r_id[147,Intercept]"
[121] "r_id[148,Intercept]" "r_id[149,Intercept]" "r_id[150,Intercept]"
[124] "r_id[151,Intercept]" "r_id[152,Intercept]" "r_id[153,Intercept]"
[127] "r_id[154,Intercept]" "r_id[155,Intercept]" "r_id[157,Intercept]"
[130] "r_id[158,Intercept]" "r_id[160,Intercept]" "r_id[162,Intercept]"
[133] "r_id[163,Intercept]" "r_id[164,Intercept]" "r_id[165,Intercept]"
[136] "r_id[166,Intercept]" "r_id[167,Intercept]" "r_id[168,Intercept]"
[139] "r_id[169,Intercept]" "r_id[170,Intercept]" "r_id[171,Intercept]"
[142] "r_id[173,Intercept]" "r_id[174,Intercept]" "r_id[176,Intercept]"
[145] "r_id[177,Intercept]" "r_id[178,Intercept]" "r_id[181,Intercept]"
[148] "r_id[182,Intercept]" "r_id[183,Intercept]" "r_id[184,Intercept]"
[151] "r_id[185,Intercept]" "r_id[186,Intercept]" "r_id[188,Intercept]"
[154] "r_id[189,Intercept]" "r_id[190,Intercept]" "r_id[191,Intercept]"
[157] "r_id[194,Intercept]" "r_id[195,Intercept]" "r_id[196,Intercept]"
[160] "r_id[197,Intercept]" "r_id[198,Intercept]" "r_id[199,Intercept]"
[163] "r_id[201,Intercept]" "r_id[202,Intercept]" "r_id[203,Intercept]"
[166] "r_id[204,Intercept]" "r_id[205,Intercept]" "r_id[206,Intercept]"
[169] "r_id[207,Intercept]" "r_id[208,Intercept]" "r_id[209,Intercept]"
[172] "r_id[211,Intercept]" "r_id[212,Intercept]" "r_id[213,Intercept]"
[175] "r_id[214,Intercept]" "r_id[216,Intercept]" "r_id[217,Intercept]"
[178] "r_id[219,Intercept]" "r_id[220,Intercept]" "r_id[221,Intercept]"
[181] "r_id[222,Intercept]" "r_id[223,Intercept]" "r_id[224,Intercept]"
[184] "r_id[225,Intercept]" "r_id[226,Intercept]" "r_id[227,Intercept]"
[187] "r_id[228,Intercept]" "r_id[229,Intercept]" "r_id[230,Intercept]"
[190] "r_id[231,Intercept]" "r_id[232,Intercept]" "r_id[233,Intercept]"
[193] "r_id[234,Intercept]" "r_id[236,Intercept]" "r_id[237,Intercept]"
[196] "r_id[238,Intercept]" "r_id[239,Intercept]" "lprior"             
[199] "lp__"                "accept_stat__"       "stepsize__"         
[202] "treedepth__"         "n_leapfrog__"        "divergent__"        
[205] "energy__"
Code 
set.seed(9)
sample_id = sample(unique(d$id), replace=F, size = 20)
distinct(d, id, female, day) %>% 
  filter(id %in% sample_id) %>% 
  add_epred_draws(m4) %>% 
  mean_qi() %>% 
  ggplot( aes(x=day, y=.epred, group = as.factor(id), color=as.factor(female))) +
  geom_line() +
  scale_color_manual(values=c("#1c5253" , "#e07a5f")) +
  labs(x="day",y="PA.std", color = "female") +
  facet_wrap(~female) +
  theme(legend.position = "top")

{visibility = “hidden”}

So far, we’ve only added covariates to the estimates of the intercept. But we can have covariates on our slope parameters. These are interactions (the effect of variable X on Y changes as a function of variable Z.) And they introduce additional parameters into our model.

\[\begin{align*} \text{Level 1} &\\ \text{PA.std}_{ij} &= \alpha_i + \beta_i(\text{day}_{ij}) + \epsilon_{ij} \\ \text{Level 2} &\\ \alpha_j &= \gamma_{00} + \gamma_{00}\text{female}_j + U_{0j} \\ \beta_j &= \gamma_{10} + \gamma_{10}\text{female}_j + U_{1j} \\ \end{align*}\]

m5 <- brm(
  data=d,
  family=gaussian,
    PA.std ~ 1 + day*female + (1 + day | id), 
  prior = c( prior(normal(0, 1), class=Intercept),
             prior(normal(0, 1), class=b), 
             prior(exponential(1), class=sd),
             prior(exponential(1), class=sigma)),
  iter=5000, warmup=1000, chains=4, cores=4, seed=9,
  file=here("files/models/m71.5")
)

Note the hyperparameter, cor(Intercept, day)

m5
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: PA.std ~ 1 + day * female + (1 + day | id) 
   Data: d (Number of observations: 13033) 
  Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 1;
         total post-warmup draws = 16000

Multilevel Hyperparameters:
~id (Number of levels: 193) 
                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)          0.79      0.04     0.72     0.88 1.00     2408     5140
sd(day)                0.01      0.00     0.01     0.01 1.00     4026     6741
cor(Intercept,day)    -0.29      0.07    -0.43    -0.15 1.00     2667     5616

Regression Coefficients:
           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     -0.16      0.09    -0.32     0.02 1.00     1290     2757
day           -0.00      0.00    -0.01    -0.00 1.00     2518     4877
female         0.49      0.11     0.27     0.71 1.00     1228     2512
day:female     0.00      0.00    -0.00     0.00 1.00     2271     5046

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.56      0.00     0.56     0.57 1.00    35997    10764

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).